// sieve.cc
// Program for sieving elements for the the General Number
// Field Sieve (GNFS).
// It is based on the NTL library, which is based on GnuMP

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <string>
#include <NTL/ZZ.h>
#include <NTL/RR.h>
#include <NTL/ZZX.h>
#include <NTL/ZZXFactoring.h>
#include <NTL/ZZ_pXFactoring.h>
#include "nfsinput.h"

// Do trial division to make sure we got a real relation!
bool checkRelation(int a,int b,char* RFB,char* AFB, nfsPoly poly) {
	NTL::ZZ p,r,pnew;
	long rfb_Bound = LineCount(RFB)/2;
	long afb_Bound = LineCount(AFB)/2; 
	
	NTL::ZZ norm = NTL::to_ZZ(a)+NTL::to_ZZ(b)*poly.m;
	if(norm < 0) {
		norm = norm * NTL::to_ZZ(-1);
	}
	
	std::ifstream rfb(RFB);
	int h = 1;
	rfb >> p;
	while(norm != 1) {
		if (norm % p == 0) {
			norm = norm/p;
		}
		else {
			if(h==rfb_Bound) {
				return false;
			}
			else {
				rfb >> p;
				rfb >> p;
				h++;
			}
	
		}
	}
	rfb.close();

	norm=NormF(poly.f,NTL::to_ZZ(a),NTL::to_ZZ(b));
	if (norm < 0) {
		norm = norm * -1;
	}
	std::ifstream afb(AFB);
	h = 1;
	afb >> p;
	afb >> r;
	while(norm != 1) {
		if (norm % p == 0) {
			while ( (a+b*r)% p != 0) {
				afb >> pnew;
				afb >> r;
				h++;
				if(pnew != p) {
					return false;
				}
				p=pnew;
			}
			norm = norm/p;
		}
		else {
			if(h==afb_Bound) {
				return false;
			}
			afb >> p;
			afb >> r;
			h++;
		}
	}
	afb.close();
	
	return true;
}

// Perform linesieving with logarithms
NTL::ZZ LineSieve(long interval,char* RFB, char* AFB, char* relations,char* LINE,long QCBnum ,nfsPoly poly, NTL::ZZ rfb_prime_bound) {
	NTL::ZZ found = NTL::to_ZZ(0);
	NTL::ZZ m = poly.m;
	NTL::ZZ p,a_tmp,tmp,r,m_modp,loga;
	int a_tmp2,p2,l;
	int linewidth=2*interval+1;
	int counter;
	NTL::ZZ zero = NTL::to_ZZ(0);
	NTL::ZZ one = NTL::to_ZZ(1);
	long rfb_Bound = LineCount(RFB)/2;
	long afb_Bound = LineCount(AFB)/2; 
	NTL::ZZ numrelations=NTL::to_ZZ(rfb_Bound)+NTL::to_ZZ(afb_Bound)+NTL::to_ZZ(QCBnum) +2;
	
	//initialize line segment
	// sieve over [-interval, interval] for each b
	NTL::ZZ bzz;
	int b = 1;
	double* alog = new double[linewidth];
	double logp;
	// fuzzfactor when doing log arithmetic
	double fuzzfactor = NTL::log(rfb_prime_bound);	
	double skip_entry= (-1777.1811);

	std::cout << "fuzzfactor: " << fuzzfactor << std::endl;
    
	
	//check to see if some values have been sieved
	std::ifstream linefile;
	linefile.open(LINE,std::ios::in);
	if( linefile.is_open() )
	{
		std::cout<<"Appending to earlier sieve!"<<std::endl;
		std::ifstream linein(LINE);
		linein >> bzz;
		linein >> found;
		linein.close();
		b=NTL::to_long(bzz)+1;
	}
	linefile.close();
	
	while(found<numrelations) {
		// initialize array with a+bm
		// remember: a=i-interval
		tmp = NTL::to_ZZ(b)*m;
		for(int i = 0; i < linewidth; i++ ){
			
			loga=NTL::to_ZZ(i-interval) + tmp;
			if(loga<zero) {
				loga= (-loga);
			}
			if(loga==zero || NTL::GCD(i-interval,b)!=1) {
				alog[i]=skip_entry;
			}
			else {
				alog[i]=fuzzfactor-NTL::log(loga);
			}
		}
		
		//*******************************
		// Test for RFB smoothness
		//*******************************
		l=0;
		std::ifstream rfb_file(RFB);
		rfb_file >> p;
		rfb_file >> m_modp;
		logp = NTL::log(p);
		counter=1;
		while(counter<rfb_Bound) {
			// a must be on the form: a=-bm+kp
			a_tmp = -(NTL::to_ZZ(b)*m_modp);
			// get a into the interval
			tmp = -(a_tmp+NTL::to_ZZ(interval))/p;
			a_tmp+=tmp*p;
			a_tmp2 = NTL::to_long(a_tmp)+interval;
			p2 = NTL::to_long(p);
			if(a_tmp2 < 0) {
				a_tmp2+=p2;
			}
			while(a_tmp2 < linewidth){
				if(alog[a_tmp2]!=skip_entry && alog[a_tmp2]<0) {
					alog[a_tmp2] += logp;
					if(alog[a_tmp2]>=0){
						l++;
					}
				}
				a_tmp2+=p2;
				
			}
			rfb_file >> p;
			rfb_file >> m_modp;
			logp = NTL::log(p);
			counter++;
		}
		rfb_file.close();
		
		//*******************************
		// Test for AFB smoothness
		//*******************************
		for(int i = 0; i < linewidth; i++ ){
			loga=NormF(poly.f,NTL::to_ZZ(i-interval),NTL::to_ZZ(b));
			if(loga==zero || alog[i]<0 || alog[i]==skip_entry) {
				alog[i]=skip_entry;
			}
			else {
				if(loga<zero) {
					loga= (-loga);
				}
				alog[i]=fuzzfactor - NTL::log(loga);
			}
		}
    		std::ifstream afb_file(AFB);
		l=0;
		counter=1;
		for(;;) {
			afb_file >> p;
			afb_file >> r;
			logp = NTL::log(p);
			// a must be on the form: a=-br+kp
			a_tmp = -(NTL::to_ZZ(b)*r);
			// get a into the interval
			while(a_tmp < NTL::to_ZZ(-interval)) {
				a_tmp+=p;
			}
			while(a_tmp > NTL::to_ZZ(-interval+NTL::to_long(p))) {
				a_tmp-=p;
			}
			a_tmp2 = NTL::to_long(a_tmp)+interval;
			p2 = NTL::to_long(p);
			while(a_tmp2 < linewidth){
				if(alog[a_tmp2]!=skip_entry && alog[a_tmp2]<0) {
					alog[a_tmp2]+=logp;
					if(alog[a_tmp2]>=0){
						l++;
					}
				}
				a_tmp2+=p2;
			}
			counter++;
			if(counter>afb_Bound) {
				break;
			}
		}
		afb_file.close();
		
		//*******************************
		// save usable relations
		//*******************************
	
		std::ofstream outfile(relations,std::ios::app);
		for(int i = 0; i <= linewidth; i++ ){
			if(alog[i]>0) {
				if(checkRelation(i-interval,b,RFB,AFB,poly)) {
					outfile << i-interval << std::endl;
					outfile << b << std::endl;
					found++;
				}
			}
		}
		outfile.close();
		// save sieving line to .line-file...
		std::ofstream lineout(LINE,std::ios::out);
		lineout << NTL::to_ZZ(b) << std::endl;
		lineout << found << std::endl;
		lineout.close();

		std::cout << "Line \t" << b  << "\tverified relations:\t" << found << "/" << numrelations << std::endl;
	
		//*******************************
		// start new sieve line
		//*******************************
		b++;
	}
	delete alog;	
	return found;	
}

int main (int count, char** args){
	if(count==1) {
		std::cout << std::endl;
		std::cout << "Usage: ./pngfs_sieve init-file" << std::endl;
		std::cout << std::endl;
		return 1;
	}
	std::clock_t start,end;
	
	NFS_instance prob = readData(args[1]);
	std::string polyfile=args[1];
	polyfile+=".poly";
	nfsPoly poly = readPoly((char*) polyfile.c_str());

	std::string rationalfile=args[1];
	rationalfile+=".rfb";
	std::string algebraicfile=args[1];
	algebraicfile+=".afb";
	std::string quadraticfile=args[1];
	quadraticfile+=".qcb";
	std::string relationsfile=args[1];
	relationsfile+=".rel";
	std::string linefile=args[1];
	linefile+=".line";
	long numQCB = LineCount((char*) quadraticfile.c_str())/2;
	NTL::RR::SetPrecision(NTL::NumBits(prob.n));
	int interval = NTL::to_long(prob.C);
	start=std::clock();
        NTL::ZZ num = LineSieve(interval,(char*) rationalfile.c_str(), (char*) algebraicfile.c_str(), (char*) relationsfile.c_str(),(char*) linefile.c_str(),numQCB,poly,prob.rfb_bound); 
	end=std::clock();	
	std::cout << num << " relations saved in file: " << relationsfile << " (time spent: " << (double) (end - start)/CLOCKS_PER_SEC << " secs)" << std::endl;
}

